function C = optimiseUncertaintyKmeans(mix, data, vars)


[N D] = size(data.Y);

Y = data.Y;
Z = vars.Z;

m = mix.varposterior.m;
T = mix.varposterior.T;

C = zeros(D);
for k = 1:mix.ncentres
    w = Z(:, k); 

    if sum(w) > 0
        wY = Y.*(w*ones(1, D));
        sum_wY = sum(wY);
        sigma = Y'*wY - sum_wY'*m(k, :) - m(k, :)'*sum_wY + sum(w)*m(k, :)'*m(k, :);
    else
        sigma = zeros(D);
    end
    
    C = C + sum(w)*T(:, :, k) + sigma;
end

C = C./N;

